% Simulate data (Appendix H.1)

%
% Simulate sample of state and control variables
%
rng('default')              
T = 1e4;                    
deflect_sss2 = deflect_sss; 
options_.simul_length = T;  
options_.oo_ = oo_;         
deflect_sss2 = compute_deflected_series(M_,options_,deflect_sss2,var_list_,0);


% Extract variables of interest from simulation
pstar_1 = exp(deflect_sss2.y(pos.pstar_1) + [0; deflect_sss2.series(pos.pstar_1,1:end-1)']);
pstar_2 = exp(deflect_sss2.y(pos.pstar_2) + [0; deflect_sss2.series(pos.pstar_2,1:end-1)']);
pstar_3 = exp(deflect_sss2.y(pos.pstar_3) + [0; deflect_sss2.series(pos.pstar_3,1:end-1)']);
pstar_4 = exp(deflect_sss2.y(pos.pstar_4) + [0; deflect_sss2.series(pos.pstar_4,1:end-1)']);
pstar_5 = exp(deflect_sss2.y(pos.pstar_5) + [0; deflect_sss2.series(pos.pstar_5,1:end-1)']);
p_1     = exp(deflect_sss2.y(pos.p_1)     + [0; deflect_sss2.series(pos.p_1,1:end-1)']);
p_2     = exp(deflect_sss2.y(pos.p_2)     + [0; deflect_sss2.series(pos.p_2,1:end-1)']);
p_3     = exp(deflect_sss2.y(pos.p_3)     + [0; deflect_sss2.series(pos.p_3,1:end-1)']);
p_4     = exp(deflect_sss2.y(pos.p_4)     + [0; deflect_sss2.series(pos.p_4,1:end-1)']);
p_5     = exp(deflect_sss2.y(pos.p_5)     + [0; deflect_sss2.series(pos.p_5,1:end-1)']);
s_1     = exp(deflect_sss2.y(pos.s_1)     + [0; deflect_sss2.series(pos.s_1,1:end-1)']);
s_2     = exp(deflect_sss2.y(pos.s_2)     + [0; deflect_sss2.series(pos.s_2,1:end-1)']);
s_3     = exp(deflect_sss2.y(pos.s_3)     + [0; deflect_sss2.series(pos.s_3,1:end-1)']);
s_4     = exp(deflect_sss2.y(pos.s_4)     + [0; deflect_sss2.series(pos.s_4,1:end-1)']);
s_5     = exp(deflect_sss2.y(pos.s_5)     + [0; deflect_sss2.series(pos.s_5,1:end-1)']);
infl    = exp(deflect_sss2.y(pos.pi)      + [0; deflect_sss2.series(pos.pi,1:end-1)']);
mc      = exp(deflect_sss2.y(pos.mctilde) + [0; deflect_sss2.series(pos.mctilde,1:end-1)']);
s       = exp(deflect_sss2.y(pos.s)       + [0; deflect_sss2.series(pos.s,1:end-1)']);
tfp     = exp(deflect_sss2.y(pos.tfp)     + [0; deflect_sss2.series(pos.tfp,1:end-1)']);
i_ann   =    (deflect_sss2.y(pos.i_ann)   + [0; deflect_sss2.series(pos.i_ann,1:end-1)']);
i       =    (deflect_sss2.y(pos.i)       + [0; deflect_sss2.series(pos.i,1:end-1)']);
eps_i   =     deflect_sss2.y(pos.nu)      + [0; deflect_sss2.series(pos.nu,1:end-1)'];


% Normalize inflation to exactly equal 1 
infl    = infl/infl(1);


%
% Simulate panel of firms using the policies in the above simulation
%

% Number of firm types and number of firms
G = 5;   
N = 1e4; 

% Initialize matrices
markup_disp  = zeros(T,1);
markup_disp2 = zeros(T,1);
markup_disp3 = zeros(T,1);
p0           = NaN(N,G);

% Initialize initial prices
p0(:,1)         = pstar_1(1);
p0(:,2)         = pstar_2(1);
p0(:,3)         = pstar_3(1);
p0(:,4)         = pstar_4(1);
p0(:,5)         = pstar_5(1);
markup_disp(1)  = var(log(p0(:)));
p               = p0;

% Price adjustments
adj1 = (rand(N,T) > theta_1);
adj2 = (rand(N,T) > theta_2);
adj3 = (rand(N,T) > theta_3);
adj4 = (rand(N,T) > theta_4);
adj5 = (rand(N,T) > theta_5);


% Loop over t to compute markup dispersion for simulated prices
for t=2:T
    pold            = p;
    p(:,1)          = adj1(:,t) .* pstar_1(t) + (1-adj1(:,t)) .* pold(:,1) / infl(t); 
    p(:,2)          = adj2(:,t) .* pstar_2(t) + (1-adj2(:,t)) .* pold(:,2) / infl(t); 
    p(:,3)          = adj3(:,t) .* pstar_3(t) + (1-adj3(:,t)) .* pold(:,3) / infl(t); 
    p(:,4)          = adj4(:,t) .* pstar_4(t) + (1-adj4(:,t)) .* pold(:,4) / infl(t); 
    p(:,5)          = adj5(:,t) .* pstar_5(t) + (1-adj5(:,t)) .* pold(:,5) / infl(t); 
    markup_disp(t)  = var(log(p(:)));
end


% Estimate LP regression on simulated data
eps_mp = eps_i / std(eps_i);
Hmax   = 16;
b      = zeros(Hmax+1,1);
for h = 0:Hmax
    bb = 100*regress(markup_disp(2+h:end) -markup_disp(1:end-1-h),[eps_mp(2:end-h) ones(T-h-1,1)]);
    b(h+1,1) = bb(1);
end
irf_mp_simulated = b;


